#ifndef KERNELS_HEADER	
#define KERNELS_HEADER
#include <cmath>


#define NUM_INT 5000
#define PI 3.14159265359f
#define DO_HACK

class Kernel {
public:
	Kernel () : mKernelH(0.0), mKernelHSqr(0.0) {
#ifdef DO_HACK
		mIntVal = 0.0;
		mRecipIntVal = 0.0;
#endif
	}

	void SetKernelH( float inh ) { 
		mKernelH = inh;
		mKernelHSqr = inh*inh; 
	
		#ifdef DO_HACK
				mIntVal = mKernelH/NUM_INT;
				mRecipIntVal = 1.0f/mIntVal;
				mIntSqrVal = mKernelHSqr/NUM_INT;
				mRecipIntSqrVal = 1.0f/mIntSqrVal;
		#endif
	}

	//Evaluates the kernel function
	float KernelPoly6( float rsqr ) {
#ifdef DO_HACK
		return KernelPoly6HACK(rsqr);
#else

		static float w = 315.0/(64.0f*PI*(std::pow(mKernelH, 9)));

		if (rsqr < mKernelHSqr) {
			return w*std::pow((mKernelHSqr-rsqr), 3.0f);
		}
	
		return 0.0f;
#endif
	}

	float KernelPoly6Laplac( float r ){
		static float w1 = (std::pow(mKernelH, 4.0f));
		static float w2 = 6*std::pow(mKernelH, 2.0f);

		float lengsq = r*r;

		if (lengsq > mKernelHSqr){
			return 0.0f;
		}

		return ( -945*(w1 - w2*lengsq + 5*lengsq*lengsq ) / (32*PI*std::pow(mKernelH, 9.0f)) );	
	}

	float KernelViscosityLaplac( float r) {
		static float w = 45.0/(PI*std::pow(mKernelH, 6.0f));

		if (r > mKernelH){
			return 0.0f;
		}

		return w*(mKernelH-r);
	}

	//Evaluates the deriviate
	float KernelPoly6Grad(  float r  ) {

		static float w = (32*PI*std::pow(mKernelH, 9.0f));

		float rsq = r*r;

		if (rsq > mKernelHSqr){
			return 0.0f;
		}

		return -945 * r * ((mKernelHSqr - rsq)*(mKernelHSqr-rsq))/w;
	}


	float KernelSpikeyGrad( float r ){
#ifdef DO_HACK
		return KernelSpikeyGradHACK(r);
#else


		static float w1 = -45.0f * (1.0f / (PI*std::pow(mKernelH, 6)));

		if (r > mKernelH){
			return 0.0f;
		}

		return w1 * std::pow((mKernelH - r), 2);	
#endif
	}

	float KernelSpikey( float r ) {

		if (r > mKernelH){
			return 0.0f;
		}

		return (15.0f / (PI*std::pow(mKernelH, 6))) * std::pow((mKernelH - r), 3);
	}



private:

#ifdef DO_HACK
	float KernelPoly6HACK(float rsqr){
		static float kp6tbl[NUM_INT];
		static bool hasHacked = false;

		if (!hasHacked)
		{
			static float w = 315.0/(64.0f*PI*(std::pow(mKernelH, 9)));
			for (int i = 0; i < NUM_INT; ++i) {
				float leng = mIntSqrVal * i;
				kp6tbl[i] = w*std::pow((mKernelHSqr-leng), 3.0f);
			}
			hasHacked = true;
		}

		unsigned int hackNdx = rsqr*mRecipIntSqrVal;
		if (hackNdx < 0 || hackNdx >= NUM_INT) {
			return 0.0f;
		}

		return kp6tbl[hackNdx];
	}


	float KernelSpikeyGradHACK( float r){
		static float ksgTbl[NUM_INT];
		static bool hasHacked = false;

		unsigned int hackNdx = r*mRecipIntVal;
		if (hackNdx < 0 || hackNdx >= NUM_INT){
			//Not within our support
			return 0.0f;
		}
		if (!hasHacked){
			static float w1 = -45.0f * (1.0f / (PI*std::pow(mKernelH, 6)));
			//float intVal = NUM_INT/supportRad; //interval value.
			//calculate hack table.
			for (int i = 0; i <NUM_INT; ++i){
				float leng = mIntVal * i;
				ksgTbl[i] = w1*std::pow((mKernelH - leng),2);
			}
			hasHacked = true;
		}
		return ksgTbl[ hackNdx ];
	}

	float mIntVal;
	float mIntSqrVal;
	float mRecipIntVal;
	float mRecipIntSqrVal;


#endif

	float mKernelH;
	float mKernelHSqr;
	float mHackIntVal;


};

#endif